####################################################################
##例子6-1公式解
# T=5/12; r=0.1; sigma=0.4; X=50; S=49; B=30
# d1=(log(S/X)+(r+0.5*sigma^2)*T)/(sigma*sqrt(T))
# d2=d1-sigma*sqrt(T)
# 
# d3=(log(S/B)+(r+0.5*sigma^2)*T)/(sigma*sqrt(T))
# d4=d3-sigma*sqrt(T)
# 
# d5=(log(S/B)-(r-0.5*sigma^2)*T)/(sigma*sqrt(T))
# d6=d5-sigma*sqrt(T)
# 
# d7=(log(S*X/B/B)-(r-0.5*sigma^2)*T)/(sigma*sqrt(T))
# d8=d7-sigma*sqrt(T)
# 
# ND42 = pnorm(d4)-pnorm(d2)
# BSN75 = (B/S)^(2*r/sigma/sigma-1)*(pnorm(d7)-pnorm(d5))
# ND31 = pnorm(d3)-pnorm(d1)
# BSN86 = (B/S)^(2*r/sigma/sigma+1)*(pnorm(d8)-pnorm(d6))
# V_exact = X*exp(-r*T)*(ND42-BSN75) - S*(ND31-BSN86)
# print(V_exact)

####################################################################
##例子6-1模拟解（一条路径）
# T=5/12; r=0.1; sigma=0.4; X=50; S=49; B=30
# M=100 #时间区间分成小区间个数
# dt=T/M; p1=(r-0.5*sigma^2)*dt; p2=sigma*sqrt(dt)
# z=rnorm(M) #生成每条路径的节点上的随机数
# my_stock=rep(0,M)
# my_stock[1]=S
# for (j in 1:M) my_stock[j+1]=my_stock[j]*exp(p1+p2*z[j])
# plot(0:M,my_stock,type='l')
# if (min(my_stock)<B) my_payoff=0 else
#   my_payoff=max(0,X-my_stock[M])
# my_option = exp(-r*T)*my_payoff
# print(my_stock[M]); print(my_option)

####################################################################
##例子6-1模拟解（很多路径）
# T=5/12; r=0.1; sigma=0.4; X=50; S=49; B=30
# dt=T/M; p1=(r-0.5*sigma^2)*dt; p2=sigma*sqrt(dt)
# M=100 #时间区间分成小区间个数
# 
# N=1000 #要模拟的样本路径的条数
# my_options_N_run=rep(0,N) #记录每条路径对应的期权价格
# 
# ##调用另一个函数
# source('~/Desktop/20180718/fm2022/fm2022-slides/my_option_onerun.R')
# for (k in 1:N) my_options_N_run[k]=my_option_onerun()
# V_mc=mean(my_options_N_run)
# print(V_mc)


# ####################################################################
# ##例子6-3 几何平均亚式期权公式解
# T=1; r=0.1; sigma=0.4; X=50; S=60
# d1=(log(S/X)+(r+sigma^2/6)*T/2)/(sigma*sqrt(T/3))
# d2=d1-sigma*sqrt(T/3)
# G_asian_exact=S*exp(-((r+sigma^2/6)*T/2))*pnorm(d1)-X*exp(-r*T)*pnorm(d2)
# print('几何平均亚式期权价格解析解：')
# print(G_asian_exact)

# ####################################################################
# ##例子6-3 亚式期权蒙特卡洛方法
T=1; r=0.1; sigma=0.4; X=50; S=60
M=100 #时间区间分成小区间个数
dt=T/M; p1=(r-0.5*sigma^2)*dt; p2=sigma*sqrt(dt)

N=1000 #模拟的样本路径一共有N条
J_R=rep(0,N) #用来保存每条路径的算术平均
G_R=rep(0,N) #用来保存每条路径的几何平均
for (k in 1:N) {
  z=rnorm(M) #生成每条路径的节点上的随机数
  my_stock=S*exp(c(0,cumsum(p1+p2*z))) #使用几何布朗运动的解析解
  #plot(0:M,my_stock,type='l')
  J_R[k]=mean(my_stock) #算术平均
  G_R[k]=exp(sum(log(my_stock))/M) #几何平均
}
# 
V_J_asian_sa=exp(-r*T)*(J_R-X)*(J_R>X)
V_G_asian_sa=exp(-r*T)*(G_R-X)*(G_R>X)
V_J_asian_mc=mean(V_J_asian_sa)
V_G_asian_mc=mean(V_G_asian_sa)
print('标准蒙特卡洛方法算术平均亚式期权价格：')
print(V_J_asian_mc)
print('标准蒙特卡洛方法几何平均亚式期权价格：')
print(V_G_asian_mc)

# ####################################################################
# ##例子6-3 算术平均亚式期权蒙特卡洛方法--控制变量方法

alpha=cov(V_J_asian_sa,V_G_asian_sa)/var(V_G_asian_sa)
V_J_asian_con_sa=exp(-r*T)*(V_J_asian_sa-
      alpha*(V_G_asian_sa-exp(r*T)*G_asian_exact))
V_J_asian_con_mc=mean(V_J_asian_con_sa)
print('控制变量蒙特卡洛方法算术平均亚式期权价格：')
print(V_J_asian_con_mc)

# ####################################################################
# ##例子6-3 算术平均亚式期权蒙特卡洛方法--控制变量方法--方差缩减效果

print('标准蒙特卡洛方法算术平均亚式期权价格的样本方差：')
print(var(V_J_asian_sa))
print('控制变量蒙特卡洛方法算术平均亚式期权价格的样本方差：')
print(var(V_J_asian_con_sa))
print('控制变量技术方差缩小比例：')
print(var(V_J_asian_con_sa)/var(V_J_asian_sa))
par(mfrow=c(2,1))
hist(V_J_asian_sa,breaks=20,xlim=c(0,100))
hist(V_J_asian_con_sa,breaks=20,xlim=c(0,100))
par(mfrow=(c(1,1)))


